On the double criticality of fluids adsorbed in disordered porous media 
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The phase transition of a fluid adsorbed in a heterogeneous system is studied with two 
simple lattice gas models within the framework of a mean-field theory. Despite the differ- 
ent origin of the heterogeneity (spatial variation of binding energy or fluid coordination 
numbers) , the fluid can undergo two phase transitions if the hosting system is sufficiently 
heterogeneous. It is clearly shown that such a polymorphism in the confined fluid results 
from the successive condensations in distinct spatial regions of the host. We have found 
the precise conditions at which one two phase transitions occur. The insight gained from 
the present study allows one to understand better some recent puzzling simulation results. 

I. INTRODUCTION 

It is well known that fluids exhibit a behavior quite different from that in the bulk when adsorbed in 
porous matrices (see Ref. [1] for a recent review). The main features are the lowering of the critical tem- 
perature and shrinking of the liquid-gas coexistence region. The effect of a porous medium on the phase 
diagram has recently been the subject of intensive theoretical studies, aiming at analyzing the relevant 
microscopic mechanisms. The matrix confinement and disorder are thought to be the key factors, which 
have been studied [2] by means of computer simulations [3-7] and analytical models [8-14]. In addition, 
the pore size characterization indicates [15] a significant impact of the pore blocking and connectivity 
effects. 

Despite a significant progress in our understanding of all these processes, there remain some controver- 
sial points. One of them is the existence of the second critical point which is believed to signal the wetting 
[3] (or " precondensation" [9]) transition. The second phase transition has been detected in a series of 
independent computer simulation studies [3-6] of Lennard- Jones fluids in a variety of matrix models. 
The location of the critical point has been found to be quite structure-dependent [4-6]. Moreover, two 
matrix realizations of the same porosity led to qualitatively different criticalitics [6] (either one or two 
transitions). Sarkisov and Monson [4] even questioned the very existence of the second critical point, 
showing that it can disappear after the averaging over several matrix realizations. 

Analytical models have not given a definite answer. Kierlik et al [9] studied the model of lattice gas 
with quenched impurities (LGQI). A mean- field approximation resulted in the usual criticality, while 
the inclusion of correlation effects (in the mean-spherical approximation) yielded the second transition. 
Unfortunately, it was not possible to draw a clear conclusion because for three-dimensional lattices the 
same approach gives the second transition even for the bulk fluid. A similar dilemma resulted from the 
theoretical approximations has been reported [11] for an off-lattice model studied within the Ornstein- 
Zernike approach. Recent rigorous results for LGQI model on a Bethe lattice [14] did not give any 
evidence for the second transition (at least in the range of parameters the authors have studied). 

Our main objective is to find out the criterion discriminating these two types of behavior. For this 
purpose, we analyze a lattice gas model in the mean-field approximation. This analytically tractable 
model incorporates all the key features listed above. For this reason, one may expect that such an 
approach can help to understand better the double criticality. 

II. ROLE OF DISORDER AND POROSITY 

A. Random binding lattice gas 

Let us consider a lattice of N sites with the coordination number equal to q. Each site can be occupied 
by fluid or matrix particles with the corresponding occupation numbers tj = (fluid particle absent at 
site i), ti = 1 (fluid particle present at site i), x% = (matrix particle present at site i) and Xi = 1 
(matrix particle absent at site i). If the number of matrix sites is N m , then the "nominal" number of 
the fluid-accessible sites is Nf = N — N m and the overall porosity is <E> = (1 — N m )/N — Nt/N. The 



matrix configuration {xi} is quenched. Different choices for the set {xt} correspond to different matrix 
realizations. The Hamiltonian is given by 



H = -W^2t i t j -KY / t i (l-x j ), (f) 

<ij> <ij> 

where W is the fluid-fluid interaction constant and K describes the fluid-matrix interaction. The first 
double sum is over the fluid- accessible nearest neighbors (xi = Xj = 1). The second summation runs over 
the fluid-matrix nearest neighbors. Such a version of the lattice gas with quenched impurities [9,10,12] is 

studied in rcf. [14]. In the limit of very high porosity, Nf w N and X)i=fi ~ 12iLi- Note that the term 
K^2j(l — xj) can be interpreted as a random energy parameter Ei associated with the site i. Such a 
random binding arises from the random environment of adsorption sites in the presence of the matrix. 
Therefore, the effect of matrix can be considered to produce simply the random binding through the pore 
space (sites occupied by fluid particles). In this way, we arrive at the random binding lattice gas model 
with the following Hamiltonian, 

H = -wj2 t i t i-J2 Eiti - ( 2 ) 

<ij> i 

By introducing the spin- like variables Sj = 2tj — 1, the above model reduces readily to the random field 
Ising model (RFIM) [8]. 

In order to study how the phase behavior depends on the averaging over the sample realizations, we 
do not focus on a specific matrix realization but consider an ensemble of matrix samples from which 
we can extract a probability distribution f{E{) for the binding energies. For simplicity, we choose a 
non-correlated multi-modal form, 

f(E i ) = ^2c a 5{E i -e a ), (3) 

a 

where c a is the concentration of sites with the binding energy e . In terms of fluctuation ti = 9 + (ti — 9), 
the product of occupation numbers can be rewritten as, 

tit, = 2 + 6(U - 9) + 9{t -9) + (U - 6){tj - 9), (4) 

where 9 is the mean fluid density whose expression will be given shortly later (see eq. (9)). Neglecting 
the product of fluctuations, i.e., the last term on the right hand side of eq. (4), we obtain, 

Utj = t i 6 + t j 6-0 2 . (5) 

Substituting eq. (5) into eq. (2) leads finally to the Hamiltonian of the mean field approximation (MFA). 
The corresponding grand partition function for a particular matrix configuration is given by, 

Nf 

E({Ei}) = e-^e^Sf u = e 'f^^N } /2 TJ [ x + e 0(„+E i+q we)^ ^ (6) 

{U} i=l 

where (3 = l/(kT) (T: temperature and k: Boltzmann constant). In the following, we consider attractive 
fluid- fluid interactions (i.e.,W > 0). 

All the thermodynamic quantities can be obtained after averaging over the quenched disorder. The 
pressure is given by 

kT r Nf 

P = TT / IT dEif(Ei) In [S({^})] - -qW9 2 /2 + kT V c a In 1 + e ^+^+iwe) (?) 

N f J 1 1 a 1 J 

For a particular disorder configuration, the fluid density varies from site to site. The density at site j is 
given by, 

< ^' )> ~S({^})^ t ' e 6 ~ 1 + e 0(»+E j+q we)- ( 8 ) 



Averaging over the disorder, we obtain the mean fluid density, 



= I dE if ( Ej ) < tj ( Ej ) >= E c » x + e ^ + e a+g we) ■ 



(9) 



It is easy to check that the above MFA results are thermodynamically consistent, i.e., P and 9 given in 
eqs. (7) and (9) satisfy the Gibbs-Duhem relation, 



zap 



= 9. 



(10) 



Considering a bimodal distribution with the binding energies eo an d ei = eo + A and concentrations 
Co and c\ = 1 — cq we arrive at 



= c t + cm, 



(11) 



where 



To = 



e P(ti+e +qW6) 
1 + e l3(ti+e +qWe) ' 



n 



e 0{n+e o +A+qW6) 
1 -|- e 0()J.+e o +A+qWe) 



are the local densities in the regions with different binding energies. Solving eq. (9) with respect to \i for 
the bimodal case, we obtain 



fj, = -e Q - qW9 + fcTln 



G - y/G 2 - A9R{9 - 1) 
2R(0 - 1) 



(12) 



where 



G = -9(1 + R) + R(l - c ) + c ; 



R 



In the limit of vanishing heterogeneity A = 0, R = 1 we recover the conventional mean-field result for a 
bulk fluid, 



H = -e - qW9 + kTln 
The compressibility can be found in the standard way, 



1-9 



X 



Solving x 1 = with respect to T, we can determine the spinodal curve. The result for the spinodal is 
analytical but requires a very heavy algebra. Moreover, the criticality conditions, 



^\ =0 
89 1 



\89 2 J, 



= 0, 



(13) 



lead to equations for T c and 9 C which are not analytically soluble. For this reason, we introduce some 
further approximations which allow us to study more explicitly the critical point in a few limiting cases. 

In the case of weak heterogeneity (small A or R — 1), we can expand (i in eq. (12) to the second order 
in (R — 1) and obtain, 



(j, = -e - qW9 + kT 



In 



9 
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(c -l)(R-l) + y(c ,9){R-lf 



where ^(cq, 0) = h — c (l — 9) — Cq(9 — |). From the criticality conditions (13), we get 



= 1 - 



Cq(1 - C ) 



,A/kT c _ x 



= 1A 



(14) 



(15) 



where T^ ulk = qW/ik is the bulk critical temperature. To the second order of (R-l), the result for the 
critical density, i.e., 9 C = 1/2, is identical to that for a bulk fluid. However, eq. (15) shows clearly that the 
critical temperature is always lower than that of the bulk fluid, i.e., T c < T^ ulk . The lowering of critical 



temperature of confined fluids has been observed in computer simulations [3-6]. Moreover, eq. (15) has 
only one solution for T c . Numerical solution of the complete MFA equation for T c gives also only one 
critical point for small A. 

In the case of strong heterogeneity, the occupation of regions with different binding energies occurs 
at quite different thermodynamic conditions {e.g., chemical potentials or pressures). At low chemical 
potential, the fluid fills essentially the region with lower external potential i.e., n << To (in the case of 
a large negative A). Neglecting n in eq. (11), we obtain 

(16) 

The criticality conditions (13) lead to 

9 C = co/2; T c /T^ lk = co. 

When the region of lower energy is completely filled, r = 1, the fluid density given in eq. (11) can be 
approximated as 6 = c + (1 — c )ri. This leads to 

/i= -e -qW6- A + fcTln - 

The second critical point is located at 

Oc = (1 + c )/2; T c /T c wfc = 1 - c . 

Note that Co is the concentration of the sites with lower external potential and the change A — > — A is 
equivalent to c — > 1 — c . The above approximate analysis reveals that there is only one critical point in 
the case of weak heterogeneity while two critical points in strongly heterogeneous systems. 

As an illustration, we show in Figure 1 the spinodal curves (as a function of fluid density 9) for different 
values of A. The critical point is shifted to lower temperatures and lower densities when |A| is increased. 
The picture is symmetrical with respect to = 1/2 for A > 0, this is equivalent to switching between the 
regions of higher and lower external potentials. The high-temperature critical point corresponds to the 
condensation in the the larger sublattice (i.e., that with a higher concentration). The matrix disorder 
plays only a marginal role, because essentially the same shift of the critical parameters is observed for 
fluids confined in regular geometries [16]. The second (low-temperature) critical point appears only when 
the heterogeneity exceeds some threshold, i.e., |A| > |A*|. The other familiar features, i.e., lowering 
T c (for both transitions) and the shrinking of the density gap between the coexisting phases, are also 
present. The inset displays the variation of the spinodals with the concentration Co (in the case of two 
transitions, |A| > |A*|). 

Figure 2 shows the variation of the critical densities with |A| (for A < 0) for a given c (c = 0.3). 
This figure shows clearly that the threshold, |A*|, for this case is near AqW . It is quite natural that the 
threshold, |A*|, changes with Co. The inset gives a "phase diagram" in the parameter space of |A*| and 
Co which delimits the region of one or two phase transitions. This phase diagram shows that the lowest 
threshold is for the case of Co = c\ and that |A*| becomes higher and higher when |co — Ci| is increased. 
This variation can be interpreted as follows. When |c — ci| — > 1 (either c — > 1 or c\ — > 1), the system 
becomes spatially more and more homogeneous. In this case, the heterogeneity level, sufficient for the 
double transition, can be reached only with a sufficiently high value of |A*|. 

In this subsection, we have shown that the bimodal random binding lattice gas model is capable of 
reproducing the basic features of the phase diagram for a fluid confined in a random porous medium. 
The model ignores the role of porosity but focuses on the role of the disorder in the adsorbate binding 
energies. The second critical point is shown to survive only if the heterogeneity exceeds some threshold 
(measured by |A*|). 

B. Lattice gas with quenched impurities 



-e - qW8 + In 



co -0 



(17) 



In order to take into account the matrix porosity and the variation of local fluid coordination, we 
consider now the lattice gas with quenched impurities [9,10,12,13] which is described by the following 
Hamiltonian, 



H = -W ^2 ktjXiXj - K ^ M^ 1 - x j) + tjXjil - Xi)] (18) 



<ij> <ij> 



The sums are over nearest neighbors. For this reason — x j) = Pi IS tnc number of the matrix 

nearest neighbors around the site i and J2j x j = q — Pi- MFA (see eq. (5)) leads to the following 
simplified Hamiltonian, 

W8 2 i — i' t — i' ^ — ^' 
H = —[N f q - £ ft] - qW9j2 U-(K- W9)J2 UPi- (19) 

i i i 

Note that in eq. (19) we count only the sites accessible to the fluid and Yl!i — Yld=\ denotes the summation 
over these sites. The grand partition function for a particular matrix realization is given by 

N f 

»=i 

For the lattice model considered here, the distribution of matrix nearest neighbors around a fluid site can 
be written as, 

f(Pi) = ^2c a 5 Pi ,p a , ^c = l. (21) 

a a 

Taking the average over p^, we obtain the following expression for the pressure, 

P = W E/(^)lnS({ Pi }) = ~[Q~Hi +kT^c a In {l + e ^we +{ K-we )Pa] y (22) 



P 

where p is the average nearest matrix neighbors around a site not occupied by a matrix particle, i.e., 

P = ^/(ft)^E'K = E c ^- ( 23 ) 

Pi f i a 

The average porosity of the matrix is given by <j> = 1 — p/q. From eq. (22), it can be clearly seen that P 
depends not only on p, but also on the individual p a 's. Therefore, two systems with the same average 
porosity can behave differently. 

The fluid density after the averaging over the matrix disorder is given by, 

e P[^+(q-p a )We+Kp a ] 

9 = Z^ Ca l + e P[p,+ (q-p a )W6+K Pa ] ■ ( 24 ) 
a 

The result given in eq. (24) agrees with that obtained by Sarkisov and Monson [13]. The only difference is 
that we do not work with local fluid densities. If the MFA (5) is applied locally, i.e., ttfj = ti0j+tj0i — 9i0j, 
we can recover the results of ref. [13]. It can be readily checked that the pressure given in eq. (22) and 
the density given in eq. (24) satisfy Gibbs-Duhem equation, eq. (10). 

From eq. (24), we can see why the critical temperature for a confined fluid is lower than in the bulk. 
The main reason is that the effective fluid-fluid interaction inside the matrix is weaker because of the 
decreasing fluid coordination number q — p a . Eq. (24) has quite the same structure as eq. (9). Therefore, 
we can anticipate that the results to be obtained with the model considered in this subsection will be 
very similar to those in the previous subsection. It is expected that there exist a region in the parameter 
space of p a 's (local number of matrix neighbors), c a 's (concentrations) and K (fluid-matrix interaction) 
in which the double criticality appears. To confirm this, we carried out numerical calculations again for 
a bimodal case with two coordinations po and p\ = po + A, such that the average number of matrix 
neighbors is, 

p = c po + (1 - co)pi = po + (1 - co) A. (25) 



Then, eq. (24) gives a typical bimodal result with two sublattices. As in the previous subsection, we can 
see either one or two phase transitions depending on the specific values of Cq and A. The spinodal curve 
with the double criticality is presented in Figure 3 which shows all the characteristic features we have 



seen already (lowered T c 's and narrowed coexistence gap). The spinodal curve is again symmetrical with 
respect to 9 = 1/2 when the sign of the fluid-matrix interaction, K, is changed. However, this does not 
change the critical temperature (for both transitions). Recall that the positive K describes a repulsive 
fluid-matrix interaction. In this case, it is difficult to associate one of the two phase transitions with the 
wetting phenomenon. 

No double criticality has been found from the rigorous results [14] for this model on a Bethe lattice. 
Nevertheless, it should be noted that the effective field distribution found by the authors of [14] is 
multimodal and temperature-dependent (see figures 3,4 in [14]). The distinction between the main peaks 
becomes less clear-cut with increasing temperature and site dilution. We may suppose that the range 
of parameters explored was out of the threshold, necessary for observing the double criticality. Another 
important point is that the authors of [14] have studied the case of fixed field, h = [i/2+qW/ 4, which is not 
the case for our study here. Such a constraint fixes a given thermodynamic path (pressure, temperature) 
along which the state of the fluid is varied. As a result, their field distribution changes its degree of 
asymmetry with the temperature. It is known that the degree of asymmetry is crucial in observing the 
double criticality in the RFIM language [8]. 

In this subsection, we studied a model which takes into account the presence of matrix particles. The 
explicit description of matrix particles introduces two key ingredients which play an important role in the 
fluid phase behavior, i.e., the local fluid coordination and the fluid- matrix interaction. It is remarkable 
that the mean-field result of this model is very similar to that obtained in the previous subsection for 
a model which only embodies the spatial heterogeneity. So, it is tempting and also looks plausible to 
consider the space heterogeneity as the key prerequisite of the double criticality. The real appearance 
of the double criticality is determined by the strength of the heterogeneity which is described by A* for 
both models considered above. Although the numerical values of A* differs for the two models, the main 
physics is not altered at all. The space heterogeneity can have different origins. For the model considered 
in this subsection, the double criticality is related to the coexistence of two polymorph liquid phases with 
different coordination numbers. This polymorphism is induced by the heterogeneity in the number of 
matrix nearest neighbors. 

III. ROLE OF EFFECTIVE DISTRIBUTION 

In the previous section, we have seen that a bimodal matrix distribution may lead cither to one or to 
two critical points. Now, a subtle question one may ask is whether it is possible to observe two phase 
transitions when the matrix distribution is apparently unimodal? Before answering this question, it is 
useful to discuss how the pore size distribution is usually determined. Even for rigid matrices, such a 
determination is not a straightforward task. For the models frequently used in computer simulations 
(e.g., matrix of hard spheres), the distribution of matrix particles is not related directly to a pore size 
distribution. The latter is usually determined through the tessellation [6] of the void space. It is to be 
noted that such tessellation methods do not allow for identifying closed cavities. The closed pores are 
inaccessible to the fluid in real experiments or in simulations which mimic the real fluid adsorption. In 
these cases, the matrix distribution seen by the fluid is different from that determined without excluding 
the closed pores. Recent studies [14,15] have demonstrated the importance of the pore blocking and 
connectivity effects. 

In order to analyze the role of closed cavities, one can assume that the nominal pore size distribution 
P(r) (i.e., that including the contribution of closed pores) is gaussian with the mean pore size a 2 and 
dispersion 82- Guassian distributions were indeed found with the tessellation method for the matrix 
samples used in simulations [6]. Since the closed pores are generated in the same way as the open ones, 
it can be reasonably expected that the distribution of the closed pores, p(r), is similar to the nominal 
one but characterized by a\ and 8\. If the fluid accessible pores and the closed cavities are uncorrelated, 
then a non-normalized distribution of fluid accessible pores g(r) is given by the product of the nominal 
distribution and the probability of finding an open pore ([1 — p(r))), i.e., 

g(r)=P(r) [1 - p(r)} . 

Even when P(r) and p(r) are unimodal distributions, g(r) could be bimodal. Such a situation is shown 
in Figure 4. One can see that the effective pore size distribution, g(r), becomes bimodal when there is 
a significant overlap between the peak of P(r) and that of p(r). Just as the appropriate combination of 
A and c leads to the double criticality, the appropriate combination of o\,(J2, 8\ and 82 results in an 
effective bimodal pore size distribution which can lead to the double criticality. This might be considered 
as one of possible explanations of the puzzling results found in the computer simulation [6] . 



From a more general point of view, the pore size distribution can also depend on the fluid state because 
of the pore blocking effects [15] or due to the matrix reaction to the fluid adsorption. This is evident for 
real world experiments with adsorption in high-porosity materials like aerogels [17, f 8]. Upon adsorption, 
such matrices can contract or expand in volume and thus change their pore size distribution. This 
effect is not exclusive to relatively soft gel-like matrices but also quite common for various intercalation 
compounds [f9,20]. In fact, the distribution of the pore space effectively seen by a fluid depends on both 
the matrix preparation and the fluid state. In our language, this is equivalent to say that the distribution 
is conditional, P(e) = P(e\6), which is conditioned by the fluid state characterized here by the fluid 
density, 9. For simplicity, we consider explicitly the binding energy distribution but all the conclusions 
can be easily translated into the porosity language. 

As we have seen in the previous sections, the average fluid density can be obtained from that for a 
particular disorder realization, 9(e), through the following relation, 



/ 



9 = / deP{e\6)6{e). (26) 



Now the compressibility becomes 
1 86 Id 



6» 2 d^i 6» 2 d^i 



(27) 



Here, the first term is an average over disorder and the second term takes into account the matrix reaction 
i.e., the modification of the matrix distribution induced by adsorbates. Eq. (27) can be rewritten as 



X = Xo 

where 



dP(e\9) Q/ 
de ^9— d{€) 



(28) 



n-hl*™*™' (29) 

In the case of a unimodal P (e\6) (no double criticality), a divergent xo (at some density 9 = 6* ) cor- 
responds to a phase transition resulting mainly from the intrinsic fluid properties. The matrix only 
rescales the critical parameters and the behavior of the confined fluid is somehow bulk-like. But, under 
the condition 



/ 



X is again singular (at another fluid density 9 — 9\). We observe under this condition a second transition. 
This effect is quite similar to the polymorphism in simple fluids with density-dependent pair potentials 
[21,22]. Moreover, due to the matrix reaction, we can see the criticality even if the bulk-like behavior is 
not critical (i.e., xo being finite). This offers a possibility for creating a critical state under conditions, 
different from the bulk critical point. This possibility might find interesting technological applications. 

We would like to emphasize that in fact the pore size distributions measured experimentally are the 
conditional ones. From the usual techniques, the pore size distribution is often deduced by fitting adsorp- 
tion isotherms [17]. It is well known that the results of such measurements depend on the fluid state and 
on the used probe molecules. From the theoretical point of view, such a distribution naturally appears 
when the maximum entropy approach [23] is applied to the determination of the relevant distribution 
(for e in the case considered here) from some indirect data 9 = 9(e) (the over-bar denotes the average 
over disorder). 

In this section, we described a scenario in which a nominal unimodal matrix distribution, found through 
analyzing the matrix geometry, can behave as an effectively bimodal if the accessibility of the pores to 
the fluid is taken into account. Alternatively, one can observe the double criticality even for a unimodal 
distribution, provided that it is sensitive to the fluid thermodynamic state. In this case, one critical point 
corresponds to a rescaled bulk-like transition, while the other transition is related to the matrix reaction 
to the presence of adsobates. 



IV. CONCLUSION 



The main question studied in the present work is how the heterogeneity of a hosting system is related 
to the occurrence of two critical points for a fluid adsorbed in it. For answering this question, we analyzed 
in details the role of several relevant factors. Our analysis is based on an analytically tractable lattice 
model with a bimodal distribution of binding energies or local porosities. We find that the the bimodality 
itself does not guarantee the existence of two critical points. The adsorbed fluid can undergo two phase 
transitions only when the host heterogeneity is sufficiently strong. For the models considered here, this 
is manifested through the existence of a threshold, A*, for the appearance of the double criticality. For 
two samples having (due to the preparation) the same mean binding energy e = eo + (1 — co)A or mean 
porosity <f> = 1 — [po + (1 — co)A]/q, one can observe either one or two critical points depending on the 
particular values taken by Co and A. This provides the criterion for predicting the occurrence of the double 
criticality for the lattice models considered here. Nevertheless, it does not seem to be straightforward to 
translate the criterion obtained here to the case of off-lattice models. 

Within the framework of the models studied here, the two phase transitions correspond to the conden- 
sation in the regions of different binding energy or local porosity. In the latter case, the nature of the 
double criticality is related to the coexistence of two polymorph liquid phases with different coordination 
numbers. This polymorphism is induced by the heterogeneity in the number of matrix nearest neighbors. 
The reduced fluid coordination number (compared to the bulk) leads to the well-known effects: lowering 
of T c and shrinking of the coexistence density gap. 

We argue that in order to give an appropriate analysis of simulation results in terms of a pore size 
distribution function, one should be aware that the matrix distribution seen effectively by the adsorbed 
fluid can depend on many factors, e.g., presence of closed pores, non rigidity of the matrix etc.. A 
nominal unimodal matrix distribution, found through analyzing the matrix geometry, can behave as a 
bimodal in the presence of a fluid if closed pores exist. This could be a possible explanation of recent 
puzzling simulation results [6]. Alternatively, one can observe the double criticality even in the case of 
an unimodal distribution, provided that it is sensitive to the fluid thermodynamic state. In this scenario, 
one critical point corresponds to a rescaled bulk-like transition, while the other transition is induced by 
the matrix reaction to the presence of the fluid. This double criticality mechanism is quite similar to the 
polymorphism in fluids with density-dependent interactions [21,22]. 



Finally, we can wonder also why no more than two phase transitions have been observed from simu- 
lations up to now. In principle, we can observe three critical points with the models studied here when 
a trimodal distribution is considered with co,ci,C2 and eo,ei,e2. In the light of the results obtained in 
this work, we expect to observe three phase transitions only under the conditions that Ai = eo — ei and 
A2 = ei — €2 are different enough. The same argument applies to the matrix coordination differences. But 
in this case, the above conditions can not be fulfilled easily because neither Ai (po — pi) nor A 2 (pi — P2) 
can exceed the maximal number of nearest matrix neighbors, q, which is restricted by the dimensionality 
of space and by steric effects. This makes the region of three phase transitions very narrow (even if it 
exists in the parameter phase diagram). This is probably why till now only two phase transitions were 
found. 
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FIG. 1. Spinodal curves for the random binding model. From the top to the bottom A = 0, — 1, —3, —5, q=4, 
c=0.7. The inset: c = 1, 0.9, 0.8, 0.6-from the top to the bottom, A = —7. The temperature is dimensionless 

T = T*/T* b 

FIG. 2. The critical coverage 9 C as a function of the energy difference A* = A/qW at Co = 0.7. The inset 
demonstrates how the threshold varies with the concentration Co. Horizontal dotted lines correspond to our 
analytical estimation for large A (see text). 

FIG. 3. Spinodal curves for the model with quenched impurities. The dashed curve corresponds to the bulk 
case (po = 0, A = 0), the symbols - K/qW = 10, c = 0.4, p = 2, q = 4, A = -1. 

FIG. 4. The overall pore size distribution P(r), the distribution of closed cavities p(r), and the resulting 
accessible pore distribution g(r). Both P(r) and p(r) are gaussian with 82 = 2, 02 — 6 and 5i = 0.5, o\ = 5 (see 
text). For pictorial purposes p(r) (squares) is reduced by a factor of 1/4. 
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